! This is a 2D advection example using the BiMocq scheme.
!
! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2024-03-24: Initial creation.

program bimocq_adv_2d_case

  use adv_2d_test_case_mod
  use bimocq_common_mod

  implicit none

  real, allocatable :: xmap(:,:,:)
  real, allocatable :: ymap(:,:,:)
  real, allocatable :: rho0(:,:)
  real, allocatable :: rho (:,:)
  integer, parameter :: ns = 3
  integer :: velocity_interp_type = 1     ! 1: bilinear; 2: biquadratic; 3: bicubic
  integer :: density_interp_type = 1      ! 1: bilinear; 2: biquadratic; 3: bicubic
  integer i, j
  character(30), parameter :: scheme = 'bimocq_2d'

  allocate(xmap(1-ns:nx+ns,1-ns:ny+ns,2))
  allocate(ymap(1-ns:nx+ns,1-ns:ny+ns,2))
  allocate(rho0(1-ns:nx+ns,1-ns:ny+ns))
  allocate(rho (1-ns:nx+ns,1-ns:ny+ns))

  dt = 0.8
  call adv_2d_test_case_init('slotted_cylinder', ns, rho)
  call reset()
  call output(scheme, 0, ns, nx, ny, x, y, rho)

  ! Run integration.
  print *, time_step, sum(rho(1:nx,1:ny))
  do while (time_step < nt)
    call advect_map()
    call interp_rho()
    call advance_time()
    call output(scheme, time_step, ns, nx, ny, x, y, rho)
    if (mod(time_step, 50) == 0) call reset()
    print *, time_step, sum(rho(1:nx,1:ny))
  end do

contains

  subroutine reset()

    do j = 1, ny
      do i = 1, nx
        xmap(i,j,old) = x(i)
        ymap(i,j,old) = y(j)
      end do
    end do
    call apply_bc(ns, nx, ny, xmap(:,:,old))
    call apply_bc(ns, nx, ny, ymap(:,:,old))
    xmap(1-ns:0    ,:,old) = xmap(1-ns:0    ,:,old) - xrange
    xmap(nx+1:nx+ns,:,old) = xmap(nx+1:nx+ns,:,old) + xrange
    ymap(:,1-ns:0    ,old) = ymap(:,1-ns:0    ,old) - yrange
    ymap(:,ny+1:ny+ns,old) = ymap(:,ny+1:ny+ns,old) + yrange
    if (time_step == 0) then
      rho0 = rho
    else
      rho0 = rho * sum(rho0(1:nx,1:ny)) / sum(rho(1:nx,1:ny))
    end if

  end subroutine reset

  subroutine interp_rho()

    real x0, y0
    integer i, j

    do j = 1, ny
      do i = 1, nx
        x0 = limit_coord(xmin, xmax, xrange, xmap(i,j,new))
        y0 = limit_coord(ymin, ymax, yrange, ymap(i,j,new))
        call interp_tracer(rho0, x0, y0, rho(i,j))
      end do
    end do
    call apply_bc(ns, nx, ny, rho)

  end subroutine interp_rho

  subroutine advect_map()

    real um, vm, xd, yd, dx, dy, xm, ym
    integer i, j, iter

    do j = 1, ny
      do i = 1, nx
        xd = x(i)
        yd = y(j)
        do iter = 1, 2
          call diff_coord(x(i), y(j), xd, yd, dx, dy)
          call add_coord(x(i), y(j), -0.5 * dx, -0.5 * dy, xm, ym)
          call interp_velocity(u, v, xm, ym, um, vm)
          call add_coord(x(i), y(j), -dt * um, -dt * vm, xd, yd)
        end do
        call interp_tracer(xmap(:,:,old), xd, yd, xmap(i,j,new))
        if (abs(xmap(i,j,new) - xmap(i,j,old)) > 0.5 * xrange) then
          xmap(i,j,new) = xmap(i,j,new) - sign(xrange, xmap(i,j,new) - xmap(i,j,old))
        end if
        call interp_tracer(ymap(:,:,old), xd, yd, ymap(i,j,new))
        if (abs(ymap(i,j,new) - ymap(i,j,old)) > 0.5 * yrange) then
          ymap(i,j,new) = ymap(i,j,new) - sign(yrange, ymap(i,j,new) - ymap(i,j,old))
        end if
      end do
    end do
    call apply_bc(ns, nx, ny, xmap(:,:,new))
    call apply_bc(ns, nx, ny, ymap(:,:,new))
    xmap(1-ns:0    ,:,new) = xmap(1-ns:0    ,:,new) - xrange
    xmap(nx+1:nx+ns,:,new) = xmap(nx+1:nx+ns,:,new) + xrange
    ymap(:,1-ns:0    ,new) = ymap(:,1-ns:0    ,new) - yrange
    ymap(:,ny+1:ny+ns,new) = ymap(:,ny+1:ny+ns,new) + yrange

  end subroutine advect_map

  subroutine diff_coord(x1, y1, x2, y2, dx, dy)

    real, intent(in) :: x1, y1, x2, y2
    real, intent(out) :: dx, dy

    real dx1, dy1, dx2, dy2

    dx = x1 - x2
    dy = y1 - y2
    dx1 = abs(dx)
    dy1 = abs(dy)
    dx2 = 1.0 - dx1
    dy2 = 1.0 - dy1
    if (dx1 > dx2) dx = -sign(dx2, dx)
    if (dy1 > dy2) dy = -sign(dy2, dy)

  end subroutine diff_coord

  subroutine add_coord(x1, y1, dx, dy, x2, y2)

    real, intent(in) :: x1, y1, dx, dy
    real, intent(out) :: x2, y2

    x2 = x1 + dx
    y2 = y1 + dy

    if (x2 < xmin) x2 = x2 + xrange
    if (x2 > xmax) x2 = x2 - xrange
    if (y2 < ymin) y2 = y2 + yrange
    if (y2 > ymax) y2 = y2 - yrange

  end subroutine add_coord

  subroutine interp_velocity(u, v, x0, y0, u0, v0)

    real, intent(in) :: u(1-ns:nx+ns,1-ns:ny+ns)
    real, intent(in) :: v(1-ns:nx+ns,1-ns:ny+ns)
    real, intent(in) :: x0, y0
    real, intent(out) :: u0, v0

    call interp2d(velocity_interp_type, ns, nx, ny, x, y, u, x0, y0, u0)
    call interp2d(velocity_interp_type, ns, nx, ny, x, y, v, x0, y0, v0)

  end subroutine interp_velocity

  subroutine interp_tracer(rho, x0, y0, rho0)

    real, intent(in) :: rho(1-ns:nx+ns,1-ns:ny+ns)
    real, intent(in) :: x0, y0
    real, intent(out) :: rho0

    call interp2d(density_interp_type, ns, nx, ny, x, y, rho, x0, y0, rho0)

  end subroutine interp_tracer

end program bimocq_adv_2d_case